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We propose an efficient procedure for numerically evolving the quantum dynamics of Delta Kicked 
Harmonic Oscillator. The method allows for longer and more accurate simulations of the system as 
well as a simple procedure for calculating the system's Floquet eigenstates and quasi-energies. The 
method is used to examine the dynamical behaviour of the system in cases where the ratio of the 
kicking frequency to the system's natural frequency is both rational and irrational. 

PACS numbers: 05.45.Pq, 03.65.-w, 05.45.Mt 

Kicked Hamiltonian systems are firmly established as prototypical models for studying chaos. Enormous effort has 
been spent examining models like the standard map and the kicked rotator. The study of kicked degenerate systems, 
those that do not obey the KAM theorem, have received relatively little attention but display fascinating dynamical 
£> . properties both classically and quantum mechanically. The Kicked Harmonic Oscillator (KHO) is an example of such 
a system. The system was originally proposed as a 2-dimensional model of kicked charges in a uniform magnetic field 
. 0. It has also been proposed as a model for electronic transport in semiconductor superlattices Q and atom optic 
ly-j ' modeling using ion-traps 0]. In particular, the quantum kicked harmonic oscillator, in the same way as the quantum 
kicked rotator may be simulated efficiently on a quantum information processor using the Quantum Fractional 
Fourier Transform 
. The Hamiltonian of the system may be written as 
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■ where the kicking potential V(q) in our system is given by 
P-i. V(q) = ncos(kq) (2) 



where k is a parameter related to ujho and the mass m. It is distinct from other widely studied systems because of 
the presence of two different frequencies in the Hamiltonian, namely the natural frequency of the Harmonic Oscillator 
luhO: and the kicking frequency ojk = 2tt/Tk- It is the ratio between these frequencies, ujho/^k = and the 

strength of the kicking potential, that determines the global behavior of the system [lj. It is appropriate to note here 
that the kicking potential we use is even. Using an odd kicking potential, say V(q) = /ism(kq), changes the global 
behaviour of the system dramatically . We do not attempt to analyse this case although the numerical techniques 
'i-H ' that we present later can be applied to any delta kicking potential. 

Understanding the quantum system analytically has proved to be problematic. While expressions have been found 
for the Floquet operator 0, these generally involve complicated infinite summations of Bessel functions and are 
too cumbersome to yield any important information about the dynamics of the system. Numerical simulations are 
therefore essential to understanding how the system evolves. However, in general, factors such as memory and 
computation time reduce the accuracy of simulations before much can be concluded about the eigensolutions or long 
time dynamical evolution of the system. We outline a procedure, based on the Fractional Fourier Transform, that 
allows for efficient numerical simulation of the system's evolution as well as an effective and straightforward approach 
to calculating the system's Floquet operator. We use this approach to confirm previous results for the system where 
R = 4, 5 as well as showing that delocalisation can occur when R is irrational if the kicking potential is strong enough. 
The Floquet operator for this system can be factorized into two parts: 

Ukho{Tk) — UsHo{T K )Ukick (3) 
where because of the delta function in Eq. (JIJ we can express the kick operator as 

Ukick = exp ( -^-ficos(kq) ) , (4) 
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which is clearly diagonal in the position basis and between kicks 

tt (t \ { -iu>HoT K {p 2 + q 2 )\ fK , 
U S ho{Tk) = exp I — I , (5) 

which is diagonal in the Fock state basis but not in the position basis or the momentum basis. Changing from the 
number basis to the position basis requires a matrix-vector multiplication (iV 2 operations) where N is the number 
of dimensions being used to approximate the infinite dimensional Hilbert space. However, it is well known that 
changing between the position and momentum basis can be done through Fourier Transforms. This fact is used to 
great advantage in the split step method In this method the exact operator Us ho , applied for the small duration 
At, can be approximated by 

TJ , / a i \ f-iuj HO Atp 2 \ f-w HO Atq 2 \ f-iu} HO Atp 2 \ 

U SHO (At) = exp ^ j exp ^ j exp ^ j (6) 

where we use the Fast Fourier Transform (FFT) 9], to perform the change of basis. Although this method is not 
exact, the error is 0(At) 3 1 and can therefore be made acceptably small by an appropriate choice of At. The number 
of operations needed to evolve the system from time t to time t + At is 0(N log N). However, using the stability 
condition for the Split-Step method ^(|; At < — , where / is the grid spacing, we roughly estimate the number of 
operations needed to accurately evolve the system over a finite time to be of 0(N 2 log N). The essential point is 
that we avoid the use of matrix-vector multiplications through the use of Fast Fourier Transforms. This technique is 
relatively efficient computationally and does not require large amounts of storage. 

It is possible however, to exploit certain symmetries in the Hamiltonian to evolve the system with far greater 
efficiency. The previous discussion regarded the Fourier Transform as a passive transformation, that is, one that 
changes the position basis to the momentum basis. However, for the harmonic oscillator, a Fourier Transform can 
also be regarded as an active transformation. In 1980, Namias defined the Fractional Fourier Transform in terms of 
the Hermite polynomials and showed that it is in fact the time evolution operator of the harmonic oscillator up to a 
phase Explicitly, taking ty(x,t) as the time-dependent wavefunction of the harmonic oscilator, he showed that 

*(x, t) = UsHo(t)V(x, 0) = e- u/2 F- t V(x, 0), (7) 

where the operator F-t is a generalisation of the ordinary Fourier transform [l2j. Indeed the latter is just the special 
case when t = ir/2. This result is important for our system because of the existence of Fast Fractional Fourier 
Transforms (FFrFT's) 13]. There are a number of alg orithms in the literature for performing the Fast Fractional 
Fourier Transform as efficiently as the regular FFT |l4|. This means that, in the position basis, the evolution of the 
wavefunction over a finite time may be written as 

*(x,T K ) = U K „o(T K )*(x,0) = e- iV (*)/* e -a'K/*F-. Tx *(x,0), (8) 

and can be performed in about 0(N log N) operations. This allows us to perform highly efficient numerical simulations 
of pure state quantum dynamics in Hilbert space dimensions of about 2 17 . This enables us to follow the diffusing 
quantum evolution for larger times and in greater detail. Simulations of such a size would be prohibitively costly in 
memory and processor time using other simulation techniques. 

As an additional advantage, the Floquet matrix itself can be calculated. This is a consequence of Eq. JSJ and the fact 
that the Fractional Fourier Transform matrix can be generated without much additional overhead |13(. This allows 
for quick and straightforward calculation of the Floquet eigenstates and eigenvalues using the standard LAPACK 
routines. 

One of the interesting properties of the system given in and @ , is the presence of a stochastic web that spans 
all of the classical phase space for certain values of the parameter R. Indeed the values of R for which this crystalline 
structure appears can be seen to be related to the tessellation of the phase plane, where the values of R = 3, 4, 6 
represent filling (tiling) the plane with triangles, squares or hexagons respectively (we ignore the degenerate cases of 
i? = 1,2). In the quantum mechanical problem these specific values of R have, for all values of the kicking strength, 
been shown to give rise to extended quantum eigenstates and therefore delocalisation of the quantum state over time 
000. 
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In the other cases where R is an integer, a quasi crystalline structure is observed in the classical ph ase p lane.In the 
quantum system there is a large suppression of the states quantum evolution with respect to energy • These 

results do not appear to contradict those of the tight binding prediction 01 > where extended states are predicted for 
all rational values of the frequency ratio. Indeed, there is a value of the kick strength, above which, the quantum 
system diffuses at a rate, less than, but comparable to that of the classical system 18], whose energy is known to 
grow linearly with time. 

In the case of irrational frequency ratios all traces of a crystalline structure disappear in the classical phase plane. 
The tight binding model predicts localisation for all values of the kick strength pjfl, and this has been supported 
numerically for low kick strengths [lfij ]. As far as we are aware there is no published numerical data for irrational 
frequency ratios for this system with large kick strength. 

To study the localisation properties of the system we first investigate the mean energy after each application of the 
Floquet operator. We calculate 

£(*) = M y + jhM, (9) 

as a function of the discrete time Tr- or the kick number n. A direct comparison may be made between the quantum 
state and a classical ensemble. Defining a — (q + ip)/y2H we have the wavefunction of a coherent state centered at 
the coordinates (q,p), 




The Wigner function of this state is a 2-D Gaussian that is centered about the point (q, p) [19| . The Husimi function 
is very similar to this and can be regarded as the convolution of the Wigner function with a 2-D Gaussian centered 
at the origin. In the numerical simulations we compare the quantum evolution of a coherent state with the classical 
evolution of an ensemble of normally distributed particles initially centered around (q,p) and possessing the same 
marginals as the corresponding Husimi function. 

When R = 4, the quantum system can be shown analytically to be identical to the well studied Kicked Harper 
Model (KHM) 6]. In this case the quantum system diffuses for all values of the kick strength fi, and for all initial 
conditions. To demonstrate this we center the initial state on a classical stable point (q,p) — (0,0), and note that 
the energy of the quantum state grows with time (or kick number n) , while the energy of the classical state remains 
bounded. From Fig. ^a), we see that the quantum system is able to tunnel through the separatrix which completely 
isolates classical trajectories. For completeness we also illustrate how both the quantum and classical systems diffuse 
when placed on an unstable point (q,p) = (0, w), on the classical system's stochastic web. It is worth noting that at 
fi = 0.5, the stochastic web is too small to allow classical diffusion. However we see that the quantum system uses 
the separatrix to significantly increase its rate of diffusion. 

In the case where R = 5, the classical phase plane shows a quasi-crystalline structure. From Fig. we can see 
that at low kicking strengths diffusion of the quantum system is somewhat suppressed. However, above a certain 
critical value of the kick strength /i, which depends on the initial position of the quantum state, this suppression no 
longer seems to exist. 

We see a similar type of global behavior when the frequency ratio is irrational even though the structure of the phase 
space patterns are quite different. Both the quantum and classical systems remain bounded up to different critical 
values of the kick strength fi (again the precise value of this parameter depends on the initial conditions). Above these 
values the the system begins to diffuse linearly. In the same way as before, there exists a range of kick strengths where 
the classical system diffuses and the quantum motion remain suppressed, see FigE{b)(c). It would appear that, in this 
region, the quantum state cannot flow along the same chaotic trajectories as the classical ensemble. This is opposite 
to the case for which R = 4, where at low kick strengths, the quantum system appears to exploit the infinitesimal 
stochastic web to aid its dispersion, see Fig^a)(b). 

As we have mentioned, using the definition of the Floquet operator given in Eq. (jSJ, we can easily find a finite 
dimensional matrix that can then be diagonalised by brute force. This method has allowed us find eigenvalues and 
eigenvectors of matrices up to 3000x3000. To go larger requires more RAM. This can then be used to examine the 
structure of these eigenstates in detail as well as examining the system's eigenvalue statistics. As an example we 
show that, in the case of irrational kick frequencies, the Floquet operator contains some eigenstates that are clearly 
localised in space and some that appear to be extended. The existence of these states highlights what we have already 



4 



found in the diffusion simulations. At low kick strength (say [i = 1), Fig. an initial state may have eigenstate 
components which are localised. This is not the case however if we examine the eigenstates for the system with a 
larger kick strength. When fi = 6 even the lower energy eigenstates appear to have some sort of extended nature Fig. 
[5] These figures do not show all of the individual eigenstate. We have zoomed in on a small section to emphasis the 
Husimi function's structural similarity to classical phase map. 

In conclusion, we have addressed the numerical problems associated with studying the quantum delta kicked har- 
monic oscillator. We have shown that the Fractional Fourier Transform may be used to infer crucial numerical 
information about the system. This includes improving the speed and accuracy with which we can evolve the system 
as well as giving an easy and effective way to analyse the system's Floquet eigenstates and quasi-energies. This im- 
proved numerical technique will assist in future studies of this systems dynamics using either a quantum information 
processor or the various physical systems mentioned in the introduction. 
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FIG. 1: Mean energy versus time for ujho/^k = 1/R = 1/4 and h = 1. (a), (b) n = 0.5 with initial state centered at (0,0) 
and (0,7r) respectively, (c), (d) fj, = 2 with initial state centered at (0,0) and (0, ir) respectively. A Hilbert space of dimension 
of 2 16 was used in the simulation. Solid (Dashed) curves represent the quantum (classical) system. 
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FIG. 2: Mean energy versus time for luho/^k = 1/R = 1/5 and h = 1. (a) /i = 1, (b) fj, = 2, (c) /i = 4, (d) /i = 6. The initial 
state was centered at (q,p) = (7.5,0). A Hilbert space of dimension of 2 16 was used in this simulation. Solid (Dashed) curves 
represent the quantum (classical) system. 
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FIG. 3: Mean energy versus time for luho/^k = 1/R = (V§ + l)/2, and ft = 1. (a) fx = 1, (b) fj, = 2, (c) /it = 4, (d) /i = 6. 
The initial state was centered at (q,p) = (it, 0).A Hilbert space of dimension of 2 17 was used for this simulation. Solid (Dashed) 
curves represent the quantum (classical) system. 
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FIG. 4: (a) and (b) Low energy eigenstates for the irrational system with kick strength fx = 1, plotted on a logarithmic scale. 
Below each state, (c) and (d), are their Husimi functions plotted in overlay with the classical phase space Poincare plot. A 
Hilbert space of dimension of 2 10 was used for these simulations. 
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FIG. 5: (a) and (b) Low energy eigenstates for the irrational system with kick strength /i = 6, plotted on a logarithmic scale. 
Below each state, (c) and (d), are their Husimi functions plotted in overlay with the classical phase space Poincare plot. A 
Hilbert space of dimension of 2 10 was used for these simulations. 



